//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : WavesComputer.cpp
//----------------------------------------------------------------------------

// upwind is not implemented correctly so upwind_ok is not defined!

#include <math.h>
#include <stdio.h>
#include <string>
#include <algorithm> // std::max and std::min
#ifndef _WINDOWS
#include "mermaid.h"
#endif
#include "defines.h"
#include "Array1D.h"
#include "Array2D.h"
#include "FloatUtils.h"
#include "IncomingWaves.h"
#include "WeakRefl.h"
#include "LowPassFilter.h"
#include "Solver.h"
#include "SolverFactory.h"
#include "WavesComputer.h"

#ifdef USE_CUDAWAVESCOMPUTER
#include "CudaRrbSolver.h"
#endif

#define COLUMNS  (unsigned)(Nx1 + 3)
#define ROWS     (unsigned)(Nx2 + 3)

#ifdef _WINDOWS
#pragma warning(disable: 4244) // conversion from 'double' to 'float', possible loss of data
#endif

#ifdef upwind_ok

#define SzzEdef() (m_upwind ? (U1E > 0.0 ? 0.0 :  U1E / m_dx1) :  U1E * m_dx1twice)
#define SzzWdef() (m_upwind ? (U1W > 0.0 ? -U1W / m_dx1 : 0.0) : -U1W * m_dx1twice)
#define SzzNdef() (m_upwind ? (U2N > 0.0 ? 0.0 :  U2N / m_dx2) :  U2N * m_dx2twice)
#define SzzSdef() (m_upwind ? (U2S > 0.0 ? -U2S / m_dx2 : 0.0) : -U2S * m_dx2twice)

#define SffEdef() (m_upwind ? (U1E > 0.0 ? 0.0 :  U1E / m_dx1) :  U1E * m_dx1twice)
#define SffWdef() (m_upwind ? (U1W > 0.0 ? -U1W / m_dx1 : 0.0) : -U1W * m_dx1twice)
#define SffNdef() (m_upwind ? (U2N > 0.0 ? 0.0 :  U2N / m_dx2) :  U2N * m_dx2twice)
#define SffSdef() (m_upwind ? (U2S > 0.0 ? -U2S / m_dx2 : 0.0) : -U2S * m_dx2twice)

#else

#define SzzEdef() ( U1E * m_dx1twice)
#define SzzWdef() (-U1W * m_dx1twice)
#define SzzNdef() ( U2N * m_dx2twice)
#define SzzSdef() (-U2S * m_dx2twice)

#define SffEdef() ( U1E * m_dx1twice)
#define SffWdef() (-U1W * m_dx1twice)
#define SffNdef() ( U2N * m_dx2twice)
#define SffSdef() (-U2S * m_dx2twice)

#endif

// Added by M.A. de Jong
#include <sys/time.h>
#define START_CPU_TIMER(time, iter) \
{ \
    gettimeofday(&start, NULL); \
    for (int t = 0; t < iter; t++) { \

#define STOP__CPU_TIMER(time, iter) \
    } \
    gettimeofday(&stop,  NULL); \
    time = (stop.tv_sec - start.tv_sec) + \
    1e-6 * (stop.tv_usec - start.tv_usec); \
    time /= iter; \
    }

static REAL tcg = 0;

/// constructor
WavesComputer::WavesComputer(int Nx1, ///< number of data points in x1 direction
                             int Nx2, ///< number of data points in x2 direction
                             const Array2D<REAL> &h0, ///< array with depth data
                             const Array2D<REAL> &U1, ///< array with speed data in x1 direction
                             const Array2D<REAL> &U2, ///< array with speed data in x2 direction
                             const Array2D<SHORT> &wet, ///< array with wet/dry data
                             Array2D<REAL>* pzeta, ///< array with waveheight data
                             Array2D<REAL>* pzetam1, ///< array with previous waveheight data
                             Array2D<REAL>* pphi, ///< array with speed potential data
                             Array2D<REAL>* pphim1, ///< array with previous speed potential data
                             Initialization init, ///< perform wave initialisation?
                             REAL time, ///< start time if not zero
                             REAL deltat, ///< last timestep if not zero start time
                             REAL pstart, ///< [s] fade-in time of pulses and incoming waves
                             const char *solvername, ///< name of solver to use
                             REAL dx1, ///< grid size in x1 direction
                             REAL dx2, ///< grid size in x2 direction
                             REAL g, ///< gravity
                             REAL tau_brk, ///< wave breaking parameter
                             int wvin_spec_count, ///< number of specifications for incoming waves
                             const Array1D<REAL> &wvin_spec_ampl, ///< incoming waves amplitudes
                             const Array1D<REAL> &wvin_spec_phase, ///< incoming waves phases
                             const Array1D<REAL> &wvin_spec_freq, ///< incoming waves frequencies
                             const Array1D<REAL> &wvin_spec_dir, ///< incoming waves directions
                             int bcN, ///< boundary condition north
                             int bcE, ///< boundary condition east
                             int bcS, ///< boundary condition south
                             int bcW, ///< boundary condition west
                             int somfeld, ///< absorb. bc type: 1=classical Sommerfeld, 2=improved
                             REAL k_wrbc, ///< weak reflection parameter
                             VariationalModel varmod, ///< variational model type
                             REAL kappa_c, ///< wave number
                             REAL psitol, ///< iterative solver tolerance
                             bool short_wave, ///< short wave equations?
                             bool upwind, ///< upwind scheme?
                             REAL eps, ///< ???
                             int M_lp) ///< low pass filter half width
    : m_Nx1(Nx1)
    , m_Nx2(Nx2)
    , m_dx1(dx1)
    , m_dx2(dx2)
    , m_dx1square(FloatUtils::square(dx1))
    , m_dx2square(FloatUtils::square(dx2))
    , m_dx1twice(1.0 / (2 * dx1))
    , m_dx2twice(1.0 / (2 * dx2))
    , m_g(g)
    , m_tau_brk(tau_brk)
    , m_bcN(bcN)
    , m_bcE(bcE)
    , m_bcS(bcS)
    , m_bcW(bcW)
    , m_short_wave(short_wave)
    , m_eps(eps)
    , m_upwind(upwind)
    , m_lp(M_lp)
    , m_h0(h0)
    , m_U1(U1)
    , m_U2(U2)
    , m_wet(wet)
    , m_zeta(*pzeta)
    , m_zetam1(*pzetam1)
    , m_phi(*pphi)
    , m_phim1(*pphim1)
    , m_pstart(pstart)
    , m_time(time)
    , m_dt_old(deltat)
    , m_dt (deltat)
    , m_pSolver(0)
    , m_pi(4.0 * atan(1.0))
    , m_dt_new(1)
    , m_zeta_mean(0.0)
    , m_sum_en_wgt(0.0)
    , m_zetam2(COLUMNS, ROWS) //
    , m_phim2(COLUMNS, ROWS) //
    , m_psi(COLUMNS, ROWS) //
    , m_psim1(COLUMNS, ROWS) //
    , m_psim2(COLUMNS, ROWS) //
    , m_dt_phi(COLUMNS, ROWS) //
    , m_dt_phim1(COLUMNS, ROWS) //
    , m_dt_zeta(COLUMNS, ROWS) //
    , m_dt_zetam1(COLUMNS, ROWS) //
    , m_matS(COLUMNS, ROWS) //
    , m_matC(COLUMNS, ROWS) //
    , m_matW(COLUMNS, ROWS) //
    , m_matE(COLUMNS, ROWS) //
    , m_matN(COLUMNS, ROWS) //
    // -----------------------------------// - Akash
    , m_matzS(COLUMNS, ROWS) //   Matrices for the Implicit variables - zeta and phi
    , m_matzC(COLUMNS, ROWS) //
    , m_matzW(COLUMNS, ROWS) //
    , m_matzN(COLUMNS, ROWS) //
    , m_matzE(COLUMNS, ROWS) //

    , m_matfS(COLUMNS, ROWS) //
    , m_matfC(COLUMNS, ROWS) //
    , m_matfW(COLUMNS, ROWS) //
    , m_matfN(COLUMNS, ROWS) //
    , m_matfE(COLUMNS, ROWS) //

    , m_matzfS(COLUMNS, ROWS) //   Matrices for Coupling of Implicit Variables
    , m_matzfC(COLUMNS, ROWS) //
    , m_matzfW(COLUMNS, ROWS) //
    , m_matzfN(COLUMNS, ROWS) //
    , m_matzfE(COLUMNS, ROWS) //

    , m_matfzS(COLUMNS, ROWS) //
    , m_matfzC(COLUMNS, ROWS) //
    , m_matfzW(COLUMNS, ROWS) //
    , m_matfzN(COLUMNS, ROWS) //
    , m_matfzE(COLUMNS, ROWS) //


    , m_matzpS(COLUMNS, ROWS) //   Matrices for Coupling of Implicit Variables
    , m_matzpC(COLUMNS, ROWS) //
    , m_matzpW(COLUMNS, ROWS) //
    , m_matzpN(COLUMNS, ROWS) //
    , m_matzpE(COLUMNS, ROWS) //

    , m_RHSzz(COLUMNS, ROWS) //
    , m_RHSff(COLUMNS, ROWS) //
    , m_RHSpp(COLUMNS, ROWS) //
    // --------------------------------------//
    , m_bpsi(COLUMNS, ROWS) //
    , m_alfa_brk(COLUMNS, ROWS) //
    , m_zeta_std(COLUMNS, ROWS) //
    , m_SffC(COLUMNS, ROWS) //
    , m_SRCp(COLUMNS, ROWS) //
    // System matrices
    , m_SRCz(COLUMNS, ROWS) //
    , m_SppC(COLUMNS, ROWS) //
    , m_SzzC(COLUMNS, ROWS) //  AM: Need to be changed now!
    , m_vD0(COLUMNS, ROWS) //
    , m_en_wgt(COLUMNS, ROWS) //
    , m_cN(COLUMNS) // cN(0:Nx1+2)
    , m_cS(COLUMNS) // cS(0:Nx1+2)
    , m_stoN(COLUMNS) // stoN(0:Nx1+2)
    , m_stoS(COLUMNS) // stoS(0:Nx1+2)
    , m_stozzN(COLUMNS) // stoN(0:Nx1+2)
    , m_stozzS(COLUMNS) // stoS(0:Nx1+2)
    , m_stozfN(COLUMNS) // stoN(0:Nx1+2)
    , m_stozfS(COLUMNS) // stoS(0:Nx1+2)
    , m_stoffN(COLUMNS) // stoN(0:Nx1+2)
    , m_stoffS(COLUMNS) // stoS(0:Nx1+2)
    , m_zeta_in_N_ext(COLUMNS) // zeta_in_N_ext(0:Nx1+2)
    , m_zeta_in_S_ext(COLUMNS) // zeta_in_S_ext(0:Nx1+2)
    , m_zeta_in_N_cen(COLUMNS) // zeta_in_N_cen(0:Nx1+2)
    , m_zeta_in_S_cen(COLUMNS) // zeta_in_S_cen(0:Nx1+2)
    , m_zeta_in_N_int(COLUMNS) // zeta_in_N_int(0:Nx1+2)
    , m_zeta_in_S_int(COLUMNS) // zeta_in_S_int(0:Nx1+2)
    , m_phi_in_N_ext(COLUMNS) // phi_in_N_ext(0:Nx1+2)
    , m_phi_in_S_ext(COLUMNS) // phi_in_S_ext(0:Nx1+2)
    , m_phi_in_N_cen(COLUMNS) // phi_in_N_cen(0:Nx1+2)
    , m_phi_in_S_cen(COLUMNS) // phi_in_S_cen(0:Nx1+2)
    , m_phi_in_N_int(COLUMNS) // phi_in_N_int(0:Nx1+2)
    , m_phi_in_S_int(COLUMNS) // phi_in_S_int(0:Nx1+2)
    , m_psi_in_N_ext(COLUMNS) // psi_in_N_ext(0:Nx1+2)
    , m_psi_in_S_ext(COLUMNS) // psi_in_S_ext(0:Nx1+2)
    , m_psi_in_N_cen(COLUMNS) // psi_in_N_cen(0:Nx1+2)
    , m_psi_in_S_cen(COLUMNS) // psi_in_S_cen(0:Nx1+2)
    , m_psi_in_N_int(COLUMNS) // psi_in_N_int(0:Nx1+2)
    , m_psi_in_S_int(COLUMNS) // psi_in_S_int(0:Nx1+2)
    , m_cE(ROWS) // cE(0:Nx2+2)
    , m_cW(ROWS) // cW(0:Nx2+2)
    , m_stoE(ROWS) // stoE(0:Nx2+2)
    , m_stoW(ROWS) // stoW(0:Nx2+2)
    , m_stozzE(COLUMNS) // stoN(0:Nx1+2)
    , m_stozzW(COLUMNS) // stoS(0:Nx1+2)
    , m_stozfE(COLUMNS) // stoN(0:Nx1+2)
    , m_stozfW(COLUMNS) // stoS(0:Nx1+2)
    , m_stoffE(COLUMNS) // stoN(0:Nx1+2)
    , m_stoffW(COLUMNS) // stoS(0:Nx1+2)
    , m_zeta_in_E_ext(ROWS) // zeta_in_E_ext(0:Nx2+2)
    , m_zeta_in_W_ext(ROWS) // zeta_in_W_ext(0:Nx2+2)
    , m_zeta_in_E_cen(ROWS) // zeta_in_E_cen(0:Nx2+2)
    , m_zeta_in_W_cen(ROWS) // zeta_in_W_cen(0:Nx2+2)
    , m_zeta_in_E_int(ROWS) // zeta_in_E_int(0:Nx2+2)
    , m_zeta_in_W_int(ROWS) // zeta_in_W_int(0:Nx2+2)
    , m_phi_in_E_ext(ROWS) // phi_in_E_ext(0:Nx2+2)
    , m_phi_in_W_ext(ROWS) // phi_in_W_ext(0:Nx2+2)
    , m_phi_in_E_cen(ROWS) // phi_in_E_cen(0:Nx2+2)
    , m_phi_in_W_cen(ROWS) // phi_in_W_cen(0:Nx2+2)
    , m_phi_in_E_int(ROWS) // phi_in_E_int(0:Nx2+2)
    , m_phi_in_W_int(ROWS) // phi_in_W_int(0:Nx2+2)
    , m_psi_in_E_ext(ROWS) // psi_in_E_ext(0:Nx2+2)
    , m_psi_in_W_ext(ROWS) // psi_in_W_ext(0:Nx2+2)
    , m_psi_in_E_cen(ROWS) // psi_in_E_cen(0:Nx2+2)
    , m_psi_in_W_cen(ROWS) // psi_in_W_cen(0:Nx2+2)
    , m_psi_in_E_int(ROWS) // psi_in_E_int(0:Nx2+2)
    , m_psi_in_W_int(ROWS) // psi_in_W_int(0:Nx2+2)
    , m_phiN(COLUMNS, 3) // phiN(0:Nx1+2,0:2)
    , m_psiN(COLUMNS, 3) // psiN(0:Nx1+2,0:2)
    , m_phiS(COLUMNS, 3) // phiS(0:Nx1+2,0:2)
    , m_psiS(COLUMNS, 3) // psiS(0:Nx1+2,0:2)
    , m_QphiN(COLUMNS, 3) // m_QphiN(0:Nx1+2,0:2)
    , m_QphiS(COLUMNS, 3) // m_QphiS(0:Nx1+2,0:2)
    , m_QpsiN(COLUMNS, 3) // m_QpsiN(0:Nx1+2,0:2)
    , m_QpsiS(COLUMNS, 3) // m_QpsiS(0:Nx1+2,0:2)
    , m_QzetaN(COLUMNS, 3) // m_QzetaN(0:Nx1+2,0:2)
    , m_QzetaS(COLUMNS, 3) // m_QzetaS(0:Nx1+2,0:2)
    , m_zetaN(COLUMNS, 3) // zetaN(0:Nx1+2,0:2)
    , m_zetaS(COLUMNS, 3) // zetaS(0:Nx1+2,0:2)
    , m_phiE(ROWS, 3) // phiE(0:Nx2+2,0:2)
    , m_psiE(ROWS, 3) // psiE(0:Nx2+2,0:2)
    , m_phiW(ROWS, 3) // phiW(0:Nx2+2,0:2)
    , m_psiW(ROWS, 3) // psiW(0:Nx2+2,0:2)
    , m_QphiE(ROWS, 3) // m_QphiE(0:Nx2+2,0:2)
    , m_QpsiE(ROWS, 3) // m_QpsiE(0:Nx2+2,0:2)
    , m_QphiW(ROWS, 3) // m_QphiW(0:Nx2+2,0:2)
    , m_QpsiW(ROWS, 3) // m_QpsiW(0:Nx2+2,0:2)
    , m_QzetaE(ROWS, 3) // m_QzetaE(0:Nx2+2,0:2)
    , m_QzetaW(ROWS, 3) // m_QzetaW(0:Nx2+2,0:2)
    , m_zetaE(ROWS, 3) // zetaE(0:Nx2+2,0:2)
    , m_zetaW(ROWS, 3) // zetaW(0:Nx2+2,0:2)
    #ifdef USE_CUDAWAVESCOMPUTER
    , m_pCudaWavesComputer(0) // Class CudaWavesComputer
    #endif
{
#ifdef USE_CUDAWAVESCOMPUTER
    m_pCudaWavesComputer = new CudaWavesComputer(m_Nx1, m_Nx2,
                                                 m_dx1, m_dx2,
                                                 m_bcN, m_bcE, m_bcS, m_bcW,
                                                 m_g);

    // if the CudaWavesComputer is used only a CUDA solver is allowed
    // (because  memory transfers between host and device are missing,
    // which makes the program much faster)
    solvername = CUDASOLVERNAME; // see "defines.h"

#endif


    initialize(solvername,
               somfeld,
               k_wrbc,
               varmod,
               kappa_c,
               psitol,
               wvin_spec_count,
               wvin_spec_ampl,
               wvin_spec_phase,
               wvin_spec_freq,
               wvin_spec_dir,
               init);

#ifdef USE_CUDAWAVESCOMPUTER
    // Copy all data to device (initialize it with constant objects):
    m_pCudaWavesComputer->initialize(&m_U1,
                                     &m_U2,
                                     &m_h0,
                                     &m_wet,
                                     &m_vD0,
                                     &m_SRCz,
                                     &m_SRCp,
                                     &m_alfa_brk,
                                     &m_zeta_std,
                                     &m_SppC,
                                     &m_stoN,
                                     &m_stoE,
                                     &m_stoS,
                                     &m_stoW);
#endif
}

/// destructor
WavesComputer::~WavesComputer()
{
#ifdef USE_CUDAWAVESCOMPUTER
    delete m_pCudaWavesComputer; // CudaWavesComputer
#endif
    delete m_pSolver;
    delete m_pIncomingWaves;
    delete m_pWrBc;
    delete m_pLpFil;
}

/// initialization
void WavesComputer::initialize(const char *solvername, ///< solver to use
                               int somfeld, ///< absorb. bc type: 1=classical Sommerfeld, 2=improved
                               REAL k_wrbc, ///< weak reflections
                               VariationalModel varmod, ///< variational model type
                               REAL kappa_c, ///< wave number
                               REAL psitol, ///< iterative solver tolerance
                               int wvin_spec_count, ///< number of wave parameters
                               const Array1D<REAL> &wvin_spec_ampl, ///< wave amplitudes
                               const Array1D<REAL> &wvin_spec_phase, ///< waves phases
                               const Array1D<REAL> &wvin_spec_freq, ///< wave frequencies
                               const Array1D<REAL> &wvin_spec_dir, ///< wave directions
                               Initialization init) ///< initialize waves?
{
    /* limit kappa * h_0 to not too small values */
    /* kappa = max( kappa, pi/2 / h0 ) ; */
    Array2D<REAL> kappa(m_zeta.dim1(), m_zeta.dim2());
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            if (m_h0[i1][i2] > 0.0)
            {
                kappa[i1][i2] = std::max(kappa_c, m_pi / 2 / m_h0[i1][i2]);
            }
            else
            {
                kappa[i1][i2] = 1.0;
            }
        }
    }

    /* wave breaking related variables are initialized */
    m_alfa_brk = 0.0;
    m_zeta_std = 0.0;

    /* construct the system matrices for the evolution of the flow equations */
    Array2D<REAL> vN0(m_zeta.dim1(), m_zeta.dim2());
    Array2D<REAL> vM0(m_zeta.dim1(), m_zeta.dim2());

    /* for the moment, we take a regular grid, constant depth and mean velocity field */
    switch (varmod)
    {
    case parabolicVM:
        // NOTE ASR: this leads to very strange results, in the original Fortran code too!
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
            {
                m_vD0[i1][i2] = 1.0 / 3.0 * m_h0[i1][i2];
                vM0[i1][i2] = 1.0 / 3.0 * m_h0[i1][i2];
                vN0[i1][i2] = 2.0 / 15.0 * FloatUtils::third(m_h0[i1][i2]);
            }
        }
        break;

    case coshypVM:
        for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
        {
            for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
            {
                if (m_h0[i1][i2] > 0.0)
                {
                    REAL vS0 = sinh(kappa[i1][i2] * m_h0[i1][i2]);
                    REAL vC0 = cosh(kappa[i1][i2] * m_h0[i1][i2]);
                    m_vD0[i1][i2] = vC0 - vS0 / (kappa[i1][i2] * m_h0[i1][i2]);
                    vM0[i1][i2] = 0.5 * kappa[i1][i2] *
                            (vS0 * vC0 - kappa[i1][i2] * m_h0[i1][i2]);
                    vN0[i1][i2] = 0.5 / kappa[i1][i2] *
                            (2.0 * kappa[i1][i2] * m_h0[i1][i2] * FloatUtils::square(vC0)
                             + kappa[i1][i2] * m_h0[i1][i2]
                             - 3.0 * vS0 * vC0);
                }
                else
                {
                    m_vD0[i1][i2] = 0.0;
                    vM0[i1][i2] = 0.0;
                    vN0[i1][i2] = 0.0;
                }
            }
        }
        break;

    default:
        ASSERT(false, "Unknown variational model");
    }

    /* initialise system matrices */
    /* row and column 0 and m_Nx1+1, resp. 0 and m_Nx2+1, are ghost points */
    // Array2D<REAL> m_matN(m_zeta.dim1(), m_zeta.dim2());
    // Array2D<REAL> m_matE(m_zeta.dim1(), m_zeta.dim2());


    m_matE = 0.0;
    m_matW = 0.0;
    m_matN = 0.0;
    m_matS = 0.0;
    m_matC = 1.0;

    m_matfC =1.0;
    m_matfW =0.0;
    m_matfS =0.0;
    m_matfE =0.0;
    m_matfN= 0.0;

    m_matzC = 1.0;
    m_matzW = 0.0;
    m_matzS = 0.0;
    m_matzE = 0.0;
    m_matzN = 0.0;

    m_matzfC = 1.0;
    m_matzfW = 0.0;
    m_matzfS = 0.0;
    m_matzfE = 0.0;
    m_matzfN = 0.0;

    m_matfzC = 1.0;
    m_matfzW = 0.0;
    m_matfzS = 0.0;
    m_matfzE = 0.0;
    m_matfzN = 0.0;

    m_matzpC = 1.0;
    m_matzpW = 0.0;
    m_matzpS = 0.0;
    m_matzpE = 0.0;
    m_matzpN = 0.0;


    m_SzzC = 1.0;
    m_SffC = 1.0;
    m_SppC = 1.0;

    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            REAL U1E = 0.5 * (m_U1[i1    ][i2    ] + m_U1[i1 + 1][i2    ]);
            REAL U1W = 0.5 * (m_U1[i1 - 1][i2    ] + m_U1[i1    ][i2    ]);
            REAL U2N = 0.5 * (m_U2[i1    ][i2    ] + m_U2[i1    ][i2 + 1]);
            REAL U2S = 0.5 * (m_U2[i1    ][i2 - 1] + m_U2[i1    ][i2    ]);

            REAL SzzE = SzzEdef();
            REAL SzzW = SzzWdef();
            REAL SzzN = SzzNdef();
            REAL SzzS = SzzSdef();
            REAL SffE = SffEdef();
            REAL SffW = SffWdef();
            REAL SffN = SffNdef();
            REAL SffS = SffSdef();

            m_SzzC[i1][i2] = -(SzzE + SzzW + SzzN + SzzS);
            m_SffC[i1][i2] = -(SffE + SffW + SffN + SffS);

            //  m_SzzC[i1][i2] =0;
            //  m_SffC[i1][i2] =0;

            REAL mult1 = 1.0 / m_dx1square;
            REAL mult2 = 1.0 / m_dx2square;
            REAL SppE = -0.5 * (vN0[i1][i2] + vN0[i1 + 1][i2 ]) * mult1;
            REAL SppW = -0.5 * (vN0[i1][i2] + vN0[i1 - 1][i2 ]) * mult1;
            REAL SppN = -0.5 * (vN0[i1][i2] + vN0[i1][ i2 + 1]) * mult2;
            REAL SppS = -0.5 * (vN0[i1][i2] + vN0[i1][ i2 - 1]) * mult2;
            m_SppC[i1][i2] = vM0[i1][i2] - (SppE + SppW + SppN + SppS);


            // phi naar zeta
            REAL h0  = m_h0[i1][i2];
            REAL h0e = mult1 * (h0 + m_h0[i1 +1 ][i2    ]);
            REAL h0w = mult1 * (h0 + m_h0[i1 -1 ][i2    ]);
            REAL h0n = mult2 * (h0 + m_h0[i1    ][i2 + 1]);
            REAL h0s = mult2 * (h0 + m_h0[i1    ][i2 - 1 ]);
            REAL SzfE = 0.5 * h0e;
            REAL SzfW = 0.5 * h0w;
            REAL SzfN = 0.5 * h0n;
            REAL SzfS = 0.5 * h0s;
            REAL SzfC = -(SzfE + SzfW + SzfN + SzfS);

            // psi naar zeta
            REAL vd0 = m_vD0[i1][i2];
            REAL SzpE = -0.25 * h0e * (vd0 + m_vD0[i1 + 1][i2    ]);
            REAL SzpW = -0.25 * h0w * (vd0 + m_vD0[i1 - 1][i2    ]);
            REAL SzpN = -0.25 * h0n * (vd0 + m_vD0[i1    ][i2 + 1]);
            REAL SzpS = -0.25 * h0s * (vd0 + m_vD0[i1    ][i2 - 1]);
            REAL SzpC = -(SzpE + SzpW + SzpN + SzpS);

            /* Correct for dry nodes */
            if (m_SppC[i1][i2] <= 0.0)
            {
                SppE = 0.0;
                SppW = 0.0;
                SppN = 0.0;
                SppS = 0.0;
                m_SppC[i1][i2] = 1.0;
            }

            /* field equations */
            m_matE[i1][i2] = SppE;
            m_matW[i1][i2] = SppW;
            m_matN[i1][i2] = SppN;
            m_matS[i1][i2] = SppS;
            m_matC[i1][i2] = m_SppC[i1][i2];

            // AM : System matrices for the Zeta and Phi variables.
            m_matzW[i1][i2] = SzzW;
            m_matzS[i1][i2] = SzzS;
            m_matzE[i1][i2] = SzzE;
            m_matzN[i1][i2] = SzzN;
            m_matzC[i1][i2] = m_SzzC[i1][i2];

            m_matfW[i1][i2] = SffW;
            m_matfS[i1][i2] = SffS;
            m_matfE[i1][i2] = SffE;
            m_matfN[i1][i2] = SffN;
            m_matfC[i1][i2] = m_SffC[i1][i2];

            // Construction of the Stencil
            // S12
            m_matzfC[i1][i2] = SzfC;
            m_matzfW[i1][i2] = SzfW;
            m_matzfS[i1][i2] = SzfS;
            m_matzfE[i1][i2] = SzfE;
            m_matzfN[i1][i2] = SzfN;

            // S21
            m_matfzC[i1][i2] = m_g;
            //   m_matfzW[j1][j2] = 0;
            //   m_matfzS[j1][j2] = 0;

            // S13 -- Used in the correction of variable Zeta and also to construct the RHS of psi equation
            m_matzpC[i1][i2] = SzpC;
            m_matzpW[i1][i2] = SzpW;
            m_matzpS[i1][i2] = SzpS;
            m_matzpE[i1][i2] = SzpE;
            m_matzpN[i1][i2] = SzpN;

        }
    }


    /* Compute the source terms for a variable current field */
    m_SRCz = 0.0;
    m_SRCp = 0.0;

    REAL SRC_flux_z;
    REAL SRC_flux_p;
    /* first the fluxes in the x1-direction */
    for (int i1 = 0; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 1; i2 <= m_Nx2; ++i2)
        {
            REAL sU1 = 0.0;

            SRC_flux_z = 0.5 * (m_h0[i1][i2] + m_h0[i1 + 1][i2])
                    * 0.5 * ((1.0 + sU1) * m_U1[i1][i2] + (1.0 - sU1) * m_U1[i1 + 1][i2]) / m_dx1;
            SRC_flux_p = 0.5 * (m_vD0[i1][i2] + m_vD0[i1 + 1][i2]) * SRC_flux_z;
            m_SRCz[i1][ i2] = m_SRCz[i1][ i2] + SRC_flux_z;
            m_SRCz[i1 + 1][i2] = m_SRCz[i1 + 1][i2] - SRC_flux_z;
            m_SRCp[i1][ i2] = m_SRCp[i1][ i2] + SRC_flux_p;
            m_SRCp[i1 + 1][i2] = m_SRCp[i1 + 1][i2] - SRC_flux_p;
        }
    }

    /* then the fluxes in the x2-direction */
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        for (int i2 = 0; i2 <= m_Nx2; ++i2)
        {
            REAL sU2 = 0.0;

            SRC_flux_z = 0.5 * (m_h0[i1][i2] + m_h0[i1][i2 + 1])
                    * 0.5 * ((1.0 + sU2) * m_U2[i1][i2] + (1.0 - sU2) * m_U2[i1][i2 + 1]) / m_dx2;
            SRC_flux_p = 0.5 * (m_vD0[i1][i2] + m_vD0[i1][i2 + 1]) * SRC_flux_z;
            m_SRCz[i1][i2 ] = m_SRCz[i1][i2 ] + SRC_flux_z;
            m_SRCz[i1][i2 + 1] = m_SRCz[i1][i2 + 1] - SRC_flux_z;
            m_SRCp[i1][i2 ] = m_SRCp[i1][i2 ] + SRC_flux_p;
            m_SRCp[i1][i2 + 1] = m_SRCp[i1][i2 + 1] - SRC_flux_p;
        }
    }

    // construct the boundary arrays for the incoming waves
    m_pIncomingWaves = new IncomingWaves(m_Nx1, m_Nx2, m_dx1, m_dx2, m_h0, m_U1, m_U2, m_wet,
                                         kappa, m_vD0, vM0, vN0, wvin_spec_count,
                                         wvin_spec_ampl, wvin_spec_phase,
                                         wvin_spec_freq, wvin_spec_dir,
                                         m_g, varmod);

    /* AM: apply boundary conditions */

    constructBoundaryStencil();


    /* Construct the Preconditioners and Solvers*/
    m_pSolver = SolverFactory::create(solvername, &m_matC, &m_matS, &m_matW, m_Nx1, m_Nx2, psitol);


    // Coupled Solver
    m_zfSolver = new RrbAkashSolver(m_matzC, m_matzS, m_matzW,
                                    m_matfC, m_matfS, m_matfW,
                                    m_matzfC, m_matzfS, m_matzfW,
                                    m_matfzC,  m_matfzS, m_matfzW,
                                    m_Nx1, m_Nx2, psitol,m_dt);


    // Store the RHS Stencil - to be used later

    storeRHS();

    /* weight factor for energy integrations */
    m_en_wgt = 1.0;
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        m_en_wgt[0] [i2] = 0.0;
        m_en_wgt[m_Nx1 + 1][i2] = 0.0;
        m_en_wgt[1] [i2] = 0.5;
        m_en_wgt[m_Nx1] [i2] = 0.5;
    }
    for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
    {
        m_en_wgt[i1][0 ] = 0.0;
        m_en_wgt[i1][m_Nx2 + 1] = 0.0;
        m_en_wgt[i1][1 ] = 0.5;
        m_en_wgt[i1][m_Nx2 ] = 0.5;
    }
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
        {
            m_sum_en_wgt += m_en_wgt[i1][i2];
        }
    }

    m_pLpFil = new LowPassFilter(m_Nx1, m_Nx2, m_lp, m_wet);
    m_pWrBc = new WeakRefl(somfeld, m_Nx1, m_Nx2);

    if (init != noInit)
    {
        if (init == spectrumInit)
        {
            m_pIncomingWaves->startup(m_time - m_dt_old, MAXDEPTH, &m_zetam1, &m_phim1);
            m_pIncomingWaves->startup(m_time, MAXDEPTH, &m_zeta, &m_phi);
        }
        startup();
    }
    else // everything is at rest
    {
        m_time = 0.0;
        m_zeta = 0.0;
        m_zetam1 = 0.0;
        m_phi = 0.0;
        m_phim1 = 0.0;

        m_zetam2 = 0.0;
        m_phim2 = 0.0;
        m_psi = 0.0;
        m_psim1 = 0.0;
        m_psim2 = 0.0;
        m_dt_zetam1 = 0.0;
        m_dt_phim1 = 0.0;
        m_dt_zeta = 0.0;
        m_dt_phi = 0.0;
    }
}

/// compute one simulation step based upon the previous one
void WavesComputer::compute(REAL dt, ///< delta T
                            const Array2D<REAL> &pressurePulse, ///< array with pulse pattern
                            double *pCgTime) ///< time used by solver
{
    ASSERT(pCgTime, "");

    // Added by M.A. de Jong

    m_dt_old = m_dt_new; // old time step
    m_dt_new = dt;
    m_time += m_dt_new; // present time

    REAL zeta_mean = 0.0;
    REAL pslow = pSlow(m_time, m_pstart);

    m_RHSzz =0.0;
    m_RHSff =0.0;
    m_RHSpp =0.0;


    { //########################################################################


        // compute the incoming wave conditions along the boundaries
        m_pIncomingWaves->compute(m_time, pslow,
                                  m_zeta_in_N_ext, m_phi_in_N_ext, m_psi_in_N_ext,
                                  m_zeta_in_N_cen, m_phi_in_N_cen, m_psi_in_N_cen,
                                  m_zeta_in_N_int, m_phi_in_N_int, m_psi_in_N_int,
                                  m_zeta_in_E_ext, m_phi_in_E_ext, m_psi_in_E_ext,
                                  m_zeta_in_E_cen, m_phi_in_E_cen, m_psi_in_E_cen,
                                  m_zeta_in_E_int, m_phi_in_E_int, m_psi_in_E_int,
                                  m_zeta_in_S_ext, m_phi_in_S_ext, m_psi_in_S_ext,
                                  m_zeta_in_S_cen, m_phi_in_S_cen, m_psi_in_S_cen,
                                  m_zeta_in_S_int, m_phi_in_S_int, m_psi_in_S_int,
                                  m_zeta_in_W_ext, m_phi_in_W_ext, m_psi_in_W_ext,
                                  m_zeta_in_W_cen, m_phi_in_W_cen, m_psi_in_W_cen,
                                  m_zeta_in_W_int, m_phi_in_W_int, m_psi_in_W_int);
        // modifies all m_* parameters (one dimensional arrays)

        // check for wave breaking if desired
        if (m_tau_brk > 0.0) // m_tau_brk is constant for all threads
        {
            waveBreaking();
            // modifies m_zeta_std and m_alfa_brk
        }


        boundaryZetaPhi();

        // modifies m_zeta and m_phi arrays

        // Compute the RHS
        ComputeRHS(m_RHSzz , m_RHSff,pslow, pressurePulse);

        // Use the Bi-CGStab method in the first pass
        m_zfSolver->solve(&m_RHSzz, &m_RHSzz, &m_zeta, &m_phi);
        std::cout << "Entered here" << std::endl;


        boundaryZetaPhi();

        ComputeRHS(m_RHSpp,pslow);

        m_pSolver->solve(&m_RHSpp, &m_psi);

        boundaryPsi();

        // Corrent Zeta
        m_zeta = MatAdd(m_zeta,scalarmult(-1*m_dt,matvec(m_matzp,m_psi)));

        // solvePsi(pslow);
        // modifies m_bpsi array, pslow is constant for all threads


        // adjustPsi();
        // modifies m_psi array



        /* short_wave, normal is: true , set to false for shallow water eq's. */
        { //--------------------------------------------------------------------
            tcg = seconds();
        } //master--------------------------------------------------------------


        if (m_short_wave)
        {
            /* use short wave corrections */
            // bspi is rechterlid, psi is oplossing


        }
        else
        {
            /* shallow-water equations */
        }


        { //--------------------------------------------------------------------
            double now = seconds();
            *pCgTime = now - tcg;
            //printf("%15.12f-%15.12f=%15.12f\n", now, tcg, *pCgTime);
        } //master--------------------------------------------------------------

        /// Correct zeta
        //  std::cout << m_dt << std::endl;



        // shift old values of psi
        //  shiftPsi();
        // modifies m_psim2 and m_psim1 arrays


        // apply Neumann boundary conditions to psi

        // modifies m_psi array


        // compute the time-derivatives of zeta and phi in the domain interior
        //  timeDerivatives(pslow, pressurePulse);
        // modifies m_dt_zeta and m_dt_phi arrays



        // time stepping of the weakly-reflective boundary conditions
        //   weaklyReflective();
        // modifies m_zetaN, m_QzetaN, m_phiN, m_QphiN, m_psiN, m_QpsiN and
        // their East, South and West counterpart N*3 dimensional arrays


        // smoothing of dt_zeta and dt_phi if desired, to prevent odd-even decoupling
        if (m_eps > 0.0) // m_eps is constant for all threads
        {
            //      smoothZetaPhi();
            //      // modifies m_dt_zeta and m_dt_phi arrays
        }


        //   leapFrog();
        // modifies m_zeta and m_phi arrays


        // apply low-pass filter if desired
        if (m_lp > 0) // m_lp is constant for all threads
        {
            //      m_pLpFil->filter(m_zeta);
            //      m_pLpFil->filter(m_phi);
            // modifies m_zeta and m_phi arrays
        }


        // compute mean wave height
        int nx1 = m_Nx1 + 1;
        int nx2 = m_Nx2 + 1;
        for (int i1 = 0; i1 <= nx1; ++i1)
        {
            for (int i2 = 0; i2 <= nx2; ++i2)
            {
                zeta_mean += m_en_wgt[i1][i2] * m_zeta[i1][i2];
            }
        } // for ---------------------------------------------------------------

    } // #######################################################################

    m_zeta_mean = zeta_mean / m_sum_en_wgt;
} // compute

/// saves wave field for continued simulation at a later time
bool WavesComputer::save(const char* filename)
{
    FILE *fp = fopen(filename, "w");
    if (!fp)
    {
        return false;
    }
    fprintf(fp, "%.7f\n", getTime());
    fprintf(fp, "%.7f\n", getDeltaTime());
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j <= m_Nx2 + 1; ++j)
        {
            fprintf(fp, "%.20f\n", m_zeta[i][j]);
            fprintf(fp, "%.20f\n", m_zetam1[i][j]);
            fprintf(fp, "%.20f\n", m_phi[i][j]);
            fprintf(fp, "%.20f\n", m_phim1[i][j]);
        }
    }
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            fprintf(fp, "%.20f\n", m_zetaN[i][j]);
            fprintf(fp, "%.20f\n", m_zetaS[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaN[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaS[i][j]);
            fprintf(fp, "%.20f\n", m_phiN[i][j]);
            fprintf(fp, "%.20f\n", m_phiS[i][j]);
            fprintf(fp, "%.20f\n", m_QphiN[i][j]);
            fprintf(fp, "%.20f\n", m_QphiS[i][j]);
            fprintf(fp, "%.20f\n", m_psiN[i][j]);
            fprintf(fp, "%.20f\n", m_psiS[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiN[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiS[i][j]);
        }
    }
    for (int j = 0; j < m_Nx2 + 1; ++j)
    {
        for (int i = 0; i < 3; ++i)
        {
            fprintf(fp, "%.20f\n", m_zetaE[i][j]);
            fprintf(fp, "%.20f\n", m_zetaW[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaE[i][j]);
            fprintf(fp, "%.20f\n", m_QzetaW[i][j]);
            fprintf(fp, "%.20f\n", m_phiE[i][j]);
            fprintf(fp, "%.20f\n", m_phiW[i][j]);
            fprintf(fp, "%.20f\n", m_QphiE[i][j]);
            fprintf(fp, "%.20f\n", m_QphiW[i][j]);
            fprintf(fp, "%.20f\n", m_psiE[i][j]);
            fprintf(fp, "%.20f\n", m_psiW[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiE[i][j]);
            fprintf(fp, "%.20f\n", m_QpsiW[i][j]);
        }
    }
    fclose(fp);
    return true;
}

/// loads wave field from saved file
bool WavesComputer::load(const char* filename)
{
    FILE *fp = fopen(filename, "r");
    if (!fp)
    {
        return false;
    }
#ifdef USEDOUBLE
    const char format[] = "%lf\n";
#else
    const char format[] = "%f\n";
#endif
    REAL time;
    REAL deltat;
    int ignore = 0;
    ignore = fscanf(fp, format, &time);
    ignore = fscanf(fp, format, &deltat);
    ASSERT(fabs(deltat - m_dt_old) < 0.0001, "deltaT in file differs from test parameter");
    m_time = time;
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j <= m_Nx2 + 1; ++j)
        {
            ignore = fscanf(fp, format, &(m_zeta[i][j]));
            ignore = fscanf(fp, format, &(m_zetam1[i][j]));
            ignore = fscanf(fp, format, &(m_phi[i][j]));
            ignore = fscanf(fp, format, &(m_phim1[i][j]));
        }
    }
    for (int i = 0; i <= m_Nx1 + 1; ++i)
    {
        for (int j = 0; j < 3; ++j)
        {
            ignore = fscanf(fp, format, &(m_zetaN[i][j]));
            ignore = fscanf(fp, format, &(m_zetaS[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaN[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaS[i][j]));
            ignore = fscanf(fp, format, &(m_phiN[i][j]));
            ignore = fscanf(fp, format, &(m_phiS[i][j]));
            ignore = fscanf(fp, format, &(m_QphiN[i][j]));
            ignore = fscanf(fp, format, &(m_QphiS[i][j]));
            ignore = fscanf(fp, format, &(m_psiN[i][j]));
            ignore = fscanf(fp, format, &(m_psiS[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiN[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiS[i][j]));
        }
    }
    for (int j = 0; j < m_Nx2 + 1; ++j)
    {
        for (int i = 0; i < 3; ++i)
        {
            ignore = fscanf(fp, format, &(m_zetaE[i][j]));
            ignore = fscanf(fp, format, &(m_zetaW[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaE[i][j]));
            ignore = fscanf(fp, format, &(m_QzetaW[i][j]));
            ignore = fscanf(fp, format, &(m_phiE[i][j]));
            ignore = fscanf(fp, format, &(m_phiW[i][j]));
            ignore = fscanf(fp, format, &(m_QphiE[i][j]));
            ignore = fscanf(fp, format, &(m_QphiW[i][j]));
            ignore = fscanf(fp, format, &(m_psiE[i][j]));
            ignore = fscanf(fp, format, &(m_psiW[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiE[i][j]));
            ignore = fscanf(fp, format, &(m_QpsiW[i][j]));
        }
    }
    fclose(fp);
    return true;
}


/// compute initial wave condition
void WavesComputer::startup()
{
    m_dt_new = m_dt_old;

    //  shift old values of zeta and phi
    m_zetam2 = m_zetam1;
    m_phim2 = m_phim1;
    m_zetam1 = m_zeta;
    m_phim1 = m_phi;

    //    for (int w = 0; w < 10; ++w)
    //    {
    //        weaklyReflective();
    //    }
    //    boundaryZetaPhi();

    // solve psi at the previous time step
    // solvePsi(1);
    // modifies m_bpsi array

    if (m_short_wave)
    {
        /* use short wave corrections */
        // bspi is rechterlid, psi is oplossing
#ifndef USE_CUDAWAVESCOMPUTER
        //   m_pSolver->solve(&m_bpsi, &m_psi);
#endif
    }
    else
    {
        /* shallow-water equations */
    }

    // shift old values of psi
    m_psim2 = m_psim1;
    m_psim1 = m_psi;

    // compute the time-derivatives of zeta and phi in the domain interior
    Array2D<REAL> dummy(m_zeta.dim1(), m_zeta.dim2());
    dummy = 0.0; // no pulses at all
    // timeDerivatives(1, dummy);
    // modifies m_dt_zeta and m_dt_phi arrays

    // compute mean wave height
    REAL zeta_mean = 0.0;
    for (int i2 = 0; i2 <= m_Nx2 + 1; ++i2)
    {
        for (int i1 = 0; i1 <= m_Nx1 + 1; ++i1)
        {
            zeta_mean += m_en_wgt[i1][i2] * m_zetam1[i1][i2];
        }
    }
    m_zeta_mean = zeta_mean / m_sum_en_wgt;

} // startup


/// returns current simulation time
REAL WavesComputer::getTime() const
{
    return m_time;
}


/// returns current time step
REAL WavesComputer::getDeltaTime() const
{
    return m_dt_new;
}


/// determines the wave breaking parameters
void WavesComputer::waveBreaking()
{
    // in incoming waves staan sinussen op de gpu
    
    // this is done inside the compute() parallelization:
    // compute wave breaking for the entire field (parallel for)

    /* local mean and standard deviation of the surface */
    /* check for wave breaking */
    REAL eps_brk = m_dt_new / m_tau_brk;
    REAL one_minus_eps_brk = 1.0 - eps_brk;
    REAL brk = eps_brk / m_tau_brk;
    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // (m_Nx1+2)*(m_Nx2+2) calculations
#ifdef openmp1
#pragma omp for //--------------------------------------------------------------
#endif
    for (int i1 = 0; i1 <= nx1; ++i1)
    {
#ifndef indices // use pointers for performance
        const REAL* ph0   = m_h0[i1];
        const SHORT* pwet = m_wet[i1];
        REAL *palfabrk    = m_alfa_brk[i1];
        REAL *pzetastd    = m_zeta_std[i1];
        REAL *pzeta       = m_zeta[i1];
        REAL *end         = pzeta + nx2 + 2;
        for (; pzeta != end; ++pzeta)
        {
            *pzetastd = sqrt(one_minus_eps_brk * FloatUtils::square(*pzetastd) +
                             eps_brk * FloatUtils::square(*pzeta - m_zeta_mean));
            *palfabrk *= one_minus_eps_brk;
            if (*pwet > 0 && (*pzetastd > 0.18 * (*ph0 + m_zeta_mean)))
            {
                *palfabrk += brk;
            }
            ++ph0;
            ++pwet;
            ++palfabrk;
            ++pzetastd;
        }
#else // use indices
        for (int i2 = 0; i2 <= nx2; ++i2)
        {
            // save some processor time by reducing address calculations
            REAL* palfabrk = &(m_alfa_brk[i1][i2]);

            m_zeta_std[i1][i2] = sqrt(one_minus_eps_brk * FloatUtils::square(m_zeta_std[i1][i2]) +
                                      eps_brk * FloatUtils::square(m_zeta[i1][i2] - m_zeta_mean));
            *palfabrk *= one_minus_eps_brk;
            if (m_wet[i1][i2] > 0 && (m_zeta_std[i1][i2] > 0.18 * (m_h0[i1][i2] + m_zeta_mean)))
            {
                *palfabrk += brk;
            }
        }
#endif
    }
    //for-------------------------------------------------------------------
}

/// applies Neumann boundary conditions to zeta and phi
void WavesComputer::boundaryZetaPhi()
{
    // this is done inside the compute() parallelization:
    // compute the boundaries (4 parallel sections)
    // (there is no overlap in the corners)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // m_Nx1 or m_Nx2 calculations per section
#ifdef openmp1
#pragma omp sections // 4
#endif
    { //====================================================================
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* north */
            if (m_bcN == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][nx2] = m_zeta[j1][m_Nx2] + m_zeta_in_N_ext[j1] - m_zeta_in_N_cen[j1];
                    m_phi [j1][nx2] = m_phi [j1][m_Nx2] + m_phi_in_N_ext [j1] - m_phi_in_N_cen[j1];
                }
            }
            else if (m_bcN == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][nx2] = m_zetaN[j1][0] + m_zeta_in_N_ext[j1];
                    m_phi [j1][nx2] = m_phiN [j1][0] + m_phi_in_N_ext [j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* east */
            if (m_bcE == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[nx1][j2] = m_zeta[m_Nx1][j2] + m_zeta_in_E_ext[j2] - m_zeta_in_E_cen[j2];
                    m_phi [nx1][j2] = m_phi [m_Nx1][j2] + m_phi_in_E_ext [j2] - m_phi_in_E_cen[j2];
                }
            }
            else if (m_bcE == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[nx1][j2] = m_zetaE[j2][0] + m_zeta_in_E_ext[j2];
                    m_phi [nx1][j2] = m_phiE [j2][0] + m_phi_in_E_ext [j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* south */
            if (m_bcS == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][0] = m_zeta[j1][1] + m_zeta_in_S_ext[j1] - m_zeta_in_S_cen[j1];
                    m_phi [j1][0] = m_phi [j1][1] + m_phi_in_S_ext [j1] - m_phi_in_S_cen[j1];
                }
            }
            else if (m_bcS == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_zeta[j1][0] = m_zetaS[j1][0] + m_zeta_in_S_ext[j1];
                    m_phi [j1][0] = m_phiS [j1][0] + m_phi_in_S_ext [j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* west */
            if (m_bcW == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[0][j2] = m_zeta[1][j2] + m_zeta_in_W_ext[j2] - m_zeta_in_W_cen[j2];
                    m_phi [0][j2] = m_phi [1][j2] + m_phi_in_W_ext [j2] - m_phi_in_W_cen[j2];
                }
            }
            else if (m_bcW == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_zeta[0][j2] = m_zetaW[j2][0] + m_zeta_in_W_ext[j2];
                    m_phi [0][j2] = m_phiW [j2][0] + m_phi_in_W_ext [j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
    } //====================================================================
}




/// applies Neumann boundary conditions to psi
void WavesComputer::boundaryPsi()
{
    // this is done inside the compute() parallelization:
    // compute the boundaries (4 parallel sections)
    // (there is no overlap in the corners)

    int nx1 = m_Nx1 + 1;
    int nx2 = m_Nx2 + 1;

    // m_Nx1 or m_Nx2 calculations per section
#ifdef openmp1
#pragma omp sections // 4
#endif
    { //====================================================================
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* north */
            if (m_bcN == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][nx2] = m_psi[j1][m_Nx2] + m_psi_in_N_ext[j1] - m_psi_in_N_cen[j1];
                }
            }
            else if (m_bcN == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][nx2] = m_psiN[j1][0] + m_psi_in_N_ext[j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* east */
            if (m_bcE == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[nx1][j2] = m_psi[m_Nx1][j2] + m_psi_in_E_ext[j2] - m_psi_in_E_cen[j2];
                }
            }
            else if (m_bcE == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[nx1][j2] = m_psiE[j2][0] + m_psi_in_E_ext[j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* south */
            if (m_bcS == 0) /* reflective Neumann boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][0] = m_psi[j1][1] + m_psi_in_S_ext[j1] - m_psi_in_S_cen[j1];
                }
            }
            else if (m_bcS == 1) /* weakly-reflective boundary */
            {
                for (int j1 = 1; j1 <= m_Nx1; ++j1)
                {
                    m_psi[j1][0] = m_psiS[j1][0] + m_psi_in_S_ext[j1];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
#ifdef openmp1
#pragma omp section
#endif
        { //----------------------------------------------------------------
            /* west */
            if (m_bcW == 0) /* reflective Neumann boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[0][j2] = m_psi[1][j2] + m_psi_in_W_ext[j2] - m_psi_in_W_cen[j2];
                }
            }
            else if (m_bcW == 1) /* weakly-reflective boundary */
            {
                for (int j2 = 1; j2 <= m_Nx2; ++j2)
                {
                    m_psi[0][j2] = m_psiW[j2][0] + m_psi_in_W_ext[j2];
                }
            }
            else
            {
                ASSERT(false, "unknown boundary condition");
            }
        } //section---------------------------------------------------------
    } //====================================================================
}




void WavesComputer::constructBoundaryStencil()
{
    /* apply boundary conditions */

    /* north */
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        if (m_bcN == 0) /* reflective Neumann boundary */
        {
            m_stoN[i1] = m_matN[i1][m_Nx2];
            m_matC[i1][m_Nx2] = m_matC[i1][m_Nx2] + m_matN[i1][m_Nx2];
            m_matN[i1][m_Nx2] = 0.0;

            m_stoffN[i1] = m_matfN[i1][m_Nx2];
            m_matfC[i1][m_Nx2] = m_matfC[i1][m_Nx2] + m_matfN[i1][m_Nx2];
            m_matfN[i1][m_Nx2] = 0.0;

            m_stozzN[i1] = m_matzN[i1][m_Nx2];
            m_matzC[i1][m_Nx2] = m_matzC[i1][m_Nx2] + m_matzN[i1][m_Nx2];
            m_matzN[i1][m_Nx2] = 0.0;


            m_stozfN[i1] = m_matzfN[i1][m_Nx2];
            m_matzfC[i1][m_Nx2] = m_matzfC[i1][m_Nx2] + m_matzfN[i1][m_Nx2];
            m_matzfN[i1][m_Nx2] = 0.0;

            m_matfzC[i1][m_Nx2] = m_matfzC[i1][m_Nx2] + m_matfzN[i1][m_Nx2];
            m_matfzN[i1][m_Nx2] = 0.0;

           // m_matzpC[i1][m_Nx2] = m_matzpC[i1][m_Nx2] + m_matzpN[i1][m_Nx2];
           // m_matzpN[i1][m_Nx2] = 0.0;



        }
        else if (m_bcN == 1) /* weakly-reflective boundary */
        {
            m_stoN[i1] = m_matN[i1][m_Nx2];
            m_matN[i1][m_Nx2] = 0.0;
            //  m_cN[i1] = std::min(sqrt(m_g * m_h0[i1][m_Nx2]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
        //  m_cN[i1] = std::min(sqrt(m_g * m_h0[i1][m_Nx2]), sqrt(m_g / k_wrbc));
    }
    // Weakly reflective B.C. Terms
    m_zetaN = 0.0; // Not stored as the
    m_phiN = 0.0;
    m_psiN = 0.0;
    m_QzetaN = 0.0;
    m_QphiN = 0.0;
    m_QpsiN = 0.0;
    //

    /* east */
    for (int i2 = 1; i2 <= m_Nx2; ++i2)
    {
        if (m_bcE == 0) /* reflective Neumann boundary */
        {
            m_stoE[i2] = m_matE[m_Nx1][i2];
            m_matC[m_Nx1][i2] = m_matC[m_Nx1][i2] + m_matE[m_Nx1][i2];
            m_matE[m_Nx1][i2] = 0.0;

            m_stoffE[i2] = m_matfE[i2][m_Nx2];
            m_matfC[m_Nx1][i2]  = m_matfC[m_Nx1][i2]  + m_matfE[m_Nx1][i2] ;
            m_matfE[m_Nx1][i2]  = 0.0;

            m_stozzE[i2] = m_matzE[i2][m_Nx2];
            m_matzC[m_Nx1][i2]  = m_matzC[m_Nx1][i2]  + m_matzE[m_Nx1][i2] ;
            m_matzE[m_Nx1][i2]  = 0.0;


            m_stozfE[i2] = m_matzfE[i2][m_Nx2];
            m_matzfC[m_Nx1][i2]  = m_matzfC[m_Nx1][i2]  + m_matzfE[m_Nx1][i2] ;
            m_matzfE[m_Nx1][i2]  = 0.0;

            m_matfzC[m_Nx1][i2]  = m_matfzC[m_Nx1][i2] + m_matfzE[m_Nx1][i2] ;
            m_matfzE[m_Nx1][i2]  = 0.0;

        //    m_matzpC[m_Nx1][i2]  = m_matzpC[m_Nx1][i2]  + m_matzpE[m_Nx1][i2] ;
         //   m_matzpE[m_Nx1][i2]  = 0.0;


        }
        else if (m_bcE == 1) /* weakly-reflective boundary */
        {
            m_stoE[i2] = m_matE[m_Nx1][i2];
            m_matE[m_Nx1][i2] = 0.0;
            //    m_cE[i2] = std::min(sqrt(m_g * m_h0[m_Nx1][i2]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    }
    m_zetaE = 0.0;
    m_phiE = 0.0;
    m_psiE = 0.0;
    m_QzetaE = 0.0;
    m_QphiE = 0.0;
    m_QpsiE = 0.0;

    /* south */
    for (int i1 = 1; i1 <= m_Nx1; ++i1)
    {
        if (m_bcS == 0) /* reflective Neumann boundary */
        {
            m_stoS[i1] = m_matS[i1][1];
            m_matC[i1][1] = m_matC[i1][1] + m_matS[i1][1];
            m_matS[i1][1] = 0.0;


            m_stoffS[i1] = m_matfS[i1][m_Nx2];
            m_matfC[i1][1] = m_matfC[i1][1] + m_matfS[i1][1];
            m_matfS[i1][m_Nx2] = 0.0;

            m_stozzS[i1] = m_matzS[i1][m_Nx2];
            m_matzC[i1][1] = m_matzC[i1][1] + m_matzS[i1][1];
            m_matzS[i1][m_Nx2] = 0.0;

            m_stozfS[i1] = m_matzfS[i1][m_Nx2];
            m_matzfC[i1][1] = m_matzfC[i1][1] + m_matzfS[i1][1];
            m_matzfS[i1][1] = 0.0;

            m_matfzC[i1][1] = m_matfzC[i1][1] + m_matfzS[i1][1];
            m_matfzS[i1][1] = 0.0;

         //   m_matzpC[i1][1] = m_matzpC[i1][1] + m_matzpS[i1][1];
         //   m_matzpS[i1][1] = 0.0;


        }
        else if (m_bcS == 1) /* weakly-reflective boundary */
        {
            m_stoS[i1] = m_matS[i1][1];
            m_matS[i1][1] = 0.0;
            //   m_cS[i1] = std::min(sqrt(m_g * m_h0[i1][1]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    }
    m_zetaS = 0.0;
    m_phiS = 0.0;
    m_psiS = 0.0;
    m_QzetaS = 0.0;
    m_QphiS = 0.0;
    m_QpsiS = 0.0;

    /* west */
    for (int i2 = 1; i2 <= m_Nx2; ++i2)
    {
        if (m_bcW == 0) /* reflective Neumann boundary */
        {
            m_stoW[i2] = m_matW[1][i2];
            m_matC[1][i2] = m_matC[1][i2] + m_matW[1][i2];
            m_matW[1][i2] = 0.0;


            m_stoffW[i2] = m_matfW[i2][m_Nx2];
            m_matfC[1][i2] = m_matfC[1][i2] + m_matfW[1][i2];
            m_matfW[1][i2] = 0.0;

            m_stozzW[i2] = m_matzW[i2][m_Nx2];
            m_matzC[1][i2] = m_matzC[1][i2] + m_matzW[1][i2];
            m_matzW[1][i2] = 0.0;

            m_stozfW[i2] = m_matzfW[i2][m_Nx2];
            m_matzfC[1][i2]= m_matzfC[1][i2] + m_matzfW[1][i2];
            m_matzfW[1][i2] = 0.0;

            m_matfzC[1][i2] = m_matfzC[1][i2] + m_matfzW[1][i2];
            m_matfzW[1][i2] = 0.0;

        //    m_matzpC[1][i2]= m_matzpC[1][i2] + m_matzpW[1][i2];
         //   m_matzpW[1][i2] = 0.0;


        }
        else if (m_bcW == 1) /* weakly-reflective boundary */
        {
            m_stoW[i2] = m_matW[1][i2];
            m_matW[1][i2] = 0.0;
            //   m_cW[i2] = std::min(sqrt(m_g * m_h0[1][i2]), sqrt(m_g / k_wrbc));
        }
        else
        {
            ASSERT(false, "unknown boundary condition");
        }
    }
    m_zetaW = 0.0;
    m_phiW = 0.0;
    m_psiW = 0.0;
    m_QzetaW = 0.0;
    m_QphiW = 0.0;
    m_QpsiW = 0.0;

    m_matzp.push_back(m_matzpC);
    m_matzp.push_back(m_matzpS);
    m_matzp.push_back(m_matzpW);
    m_matzp.push_back(m_matzpN);
    m_matzp.push_back(m_matzpE);

    m_matpf= m_matzp;
}

void WavesComputer::storeRHS()
{

    Array2D<REAL> cc = m_matzC ;
    Array2D<REAL> cs = m_matzS;
    Array2D<REAL> cw = m_matzW;
    Array2D<REAL> cn = m_matzN;
    Array2D<REAL> ce = m_matzE;


    REAL dt1 = 1/m_dt;
    REAL beta =-0.5;
    for (int ii = 0; ii <= m_Nx1 + 2; ++ii)
    {
        for (int jj = 0; jj <= m_Nx2 + 2; ++jj)
        {

            cc[ii][jj] =dt1 + cc[ii][jj]*beta;
            cs[ii][jj]*=beta;
            cw[ii][jj]*=beta;
            cn[ii][jj]*=beta;
            ce[ii][jj]*=beta;
        }
    }

    m_matzz.push_back(cc);
    m_matzz.push_back(cs);
    m_matzz.push_back(cw);
    m_matzz.push_back(cn);
    m_matzz.push_back(ce);

    // Store the matrix of S_PhiPhi
    cc = m_matfC;
    cs = m_matfS;
    cw = m_matfW;
    cn = m_matfN;
    ce = m_matfE;


    for (int ii = 0; ii <= m_Nx1 + 2; ++ii)
    {
        for (int jj = 0; jj <= m_Nx2 + 2; ++jj)
        {

            cc[ii][jj] =dt1 + cc[ii][jj]*beta;
            cs[ii][jj]*=beta;
            cw[ii][jj]*=beta;
            cn[ii][jj]*=beta;
            ce[ii][jj]*=beta;
        }
    }

    m_matff.push_back(cc);
    m_matff.push_back(cs);
    m_matff.push_back(cw);
    m_matff.push_back(cn);
    m_matff.push_back(ce);


    // Store the matrix of S_ZetaPhi
    cc = m_matzfC;
    cs = m_matzfS;
    cw = m_matzfW;
    cn = m_matzfN;
    ce = m_matzfE;

    for (int ii = 0; ii <= m_Nx1 + 2; ++ii)
    {
        for (int jj = 0; jj <= m_Nx2 + 2; ++jj)
        {

            //cc[ii][jj] +=dt1;
            cs[ii][jj]*=beta;
            cw[ii][jj]*=beta;
            cn[ii][jj]*=beta;
            ce[ii][jj]*=beta;
        }
    }



    m_matzf.push_back(cc);
    m_matzf.push_back(cs);
    m_matzf.push_back(cw);
    m_matzf.push_back(cn);
    m_matzf.push_back(ce);

    // AM: Store the Array in the Vector Form for S_zeta_phi
    cc = m_matfzC;
    cs = m_matzfS;
    cw = m_matzfW;
    cn = m_matzfN;
    ce = m_matzfE;


    for (int ii = 0; ii <= m_Nx1 + 2; ++ii)
    {
        for (int jj = 0; jj <= m_Nx2 + 2; ++jj)
        {
            cs[ii][jj]*=beta;
            cw[ii][jj]*=beta;
            cn[ii][jj]*=beta;
            ce[ii][jj]*=beta;
        }
    }

    m_matfz.push_back(cc);
    m_matfz.push_back(cs);
    m_matfz.push_back(cw);
    m_matfz.push_back(cn);
    m_matfz.push_back(ce);

}

void WavesComputer::ComputeRHS(Array2D<REAL> &RHS_ZZ , Array2D<REAL> &RHS_FF, REAL pslow, const Array2D<REAL> &pressurePulse)
{

    int nx1 = m_Nx1-1;
    int nx2 = m_Nx2-1;

    RHS_ZZ = MatAdd(matvec(m_matzz,m_zeta), matvec(m_matzf,m_phi));
    RHS_ZZ = MatSub(RHS_ZZ,scalarmult(pslow,m_SRCz));

    RHS_FF = MatAdd(matvec(m_matfz,m_zeta), matvec(m_matff,m_phi));
    RHS_FF =MatSub(RHS_FF,pressurePulse);

    if (m_bcN == 0) /* Neumann b.c. */
    {
        for (int j1 = 2; j1 <= nx1; ++j1)
        {

            RHS_ZZ[j1][m_Nx2] -= (m_stozzN[j1] * (m_zeta_in_N_ext[j1] - m_zeta_in_N_cen[j1])  +  m_stozfN[j1] * (m_phi_in_N_ext[j1] - m_phi_in_N_cen[j1]));
            RHS_FF[j1][m_Nx2] -= ( m_stoffN[j1] * (m_phi_in_N_ext[j1] - m_phi_in_N_cen[j1]));

        }
    }

    if (m_bcE == 0) /* Neumann b.c. */
    {
        for (int j2 = 2; j2 <= nx2; ++j2)
        {

            RHS_ZZ[m_Nx1][j2] -= (m_stozzE[j2] * (m_zeta_in_E_ext[j2] - m_zeta_in_E_cen[j2])  +  m_stozfE[j2] * (m_phi_in_E_ext[j2] - m_phi_in_E_cen[j2]));
            RHS_FF[m_Nx1][j2]  -= ( m_stoffE[j2] * (m_phi_in_E_ext[j2] - m_phi_in_E_cen[j2]));

        }
    }

    if (m_bcS == 0) /* Neumann b.c. */
    {
        for (int j1 = 2; j1 <= nx1; ++j1)
        {

            RHS_ZZ[j1][1] -= (m_stozzS[j1] * (m_zeta_in_S_ext[j1] - m_zeta_in_S_cen[j1])  +  m_stozfS[j1] * (m_phi_in_S_ext[j1] - m_phi_in_S_cen[j1]));
            RHS_FF[j1][1] -= ( m_stoffS[j1] * (m_phi_in_S_ext[j1] - m_phi_in_S_cen[j1]));
        }
    }

    if (m_bcW == 0) /* Neumann b.c. */
    {
        for (int j2 = 2; j2 <= nx2; ++j2)
        {

            RHS_ZZ[1][j2] -= (m_stozzW[j2] * (m_zeta_in_W_ext[j2] - m_zeta_in_W_cen[j2])  +  m_stozfW[j2] * (m_phi_in_W_ext[j2] - m_phi_in_W_cen[j2]));
            RHS_FF[1][j2]  -= ( m_stoffW[j2] * (m_phi_in_W_ext[j2] - m_phi_in_W_cen[j2]));

        }
    }

}


void WavesComputer::ComputeRHS(Array2D<float> &RHS_PP, REAL pslow)
{

    int nx1 = m_Nx1-1;
    int nx2 = m_Nx2-1;

    RHS_PP = MatSub(matvec(m_matpf,m_phi),scalarmult(pslow,m_SRCp));


    if (m_bcN == 0) /* Neumann b.c. */
    {
        for (int j1 = 2; j1 <= nx1; ++j1)
        {
           RHS_PP[j1][m_Nx2] -= ( m_stoN[j1] * (m_psi_in_N_ext[j1] - m_psi_in_N_cen[j1]));
        }
    }

    if (m_bcE == 0) /* Neumann b.c. */
    {
        for (int j2 = 2; j2 <= nx2; ++j2)
        {

            RHS_PP[m_Nx1][j2]  -= ( m_stoE[j2] * (m_psi_in_E_ext[j2] - m_psi_in_E_cen[j2]));
        }
    }

    if (m_bcS == 0) /* Neumann b.c. */
    {
        for (int j1 = 2; j1 <= nx1; ++j1)
        {

            RHS_PP[j1][1] -= ( m_stoS[j1] * (m_psi_in_S_ext[j1] - m_psi_in_S_cen[j1]));
        }
    }

    if (m_bcW == 0) /* Neumann b.c. */
    {
        for (int j2 = 2; j2 <= nx2; ++j2)
        {

            RHS_PP[1][j2]  -= ( m_stoW[j2] * (m_psi_in_W_ext[j2] - m_psi_in_W_cen[j2]));

        }
    }


}
